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TERAHERTZ  ELECTROMAGNETIC  IMAGING  OF  DIELECTRIC 

MATERIALS 


Gabriella  A.  Pinter 
Department  of  Mathematical  Sciences 
University  of  Wisconsin,  Milwaukee 


1  Introduction 

Research  efforts  were  directed  towards  the  understanding  of  high  frequency  (gigahertz  and 
higher)  ultrashort  pulse  propagation  in  dielectric  materials  for  the  purpose  of  extending 
current  imaging  capabilities  and  the  assessment  of  safety  standards.  We  formulated  a  full- 
wave  variational  Maxwell  model  that  provides  the  theoretical  and  computational  foundation 
for  this  problem.  The  model  can  incorporate  linear  and  nonlinear  polarization  dynamics,  and 
facilitates  the  comparison  of  the  resulting  phenomena.  In  particular,  we  were  interested  in 
temporal  transients,  like  the  Brillouin  precursors,  and  how  they  changed  their  characteristics 
under  different  constitutive  laws. 

The  results  described  first  concern  the  theoretical  foundations  of  this  problem,  and  show 
that  under  fairly  general  conditions  the  model  permits  a  unique  weak  solution  that  depends 
continuously  on  the  input  parameters. 

An  efficient  numerical  algorithm  was  developed  for  our  extensive  simulations  and  analysis 
based  on  different  versions  of  the  presented  model.  This  algorithm  can  be  extended  readily 
to  include  the  imaging  problem,  i.e.,  the  identification  of  electromagnetic  and  geometric 
properties  of  a  hidden  substance. 

We  also  investigated  the  interaction  of  electromagnetic  waves  with  the  dynamics  of  poly¬ 
meric/viscoelastic  materials.  As  a  first  step,  we  developed  a  molecular  based  multiscale 
model  that  correctly  describes  the  characteristics  (e.g.,  hysteresis)  exhibited  by  these  mate¬ 
rials  in  dynamical  extension  and  shear  deformations.  This  model  can  then  be  coupled  to  the 
electromagnetic  phenomena  to  provide  a  complete  description  of  the  problem.  This  research 
was  joint  work  with  H.T.  Banks  and  N.G.  Medhin  from  North  Carolina  State  University. 


2  Theoretical  results 

We  consider  a  time-domain  model  for  the  propagation  of  high-intensity  electromagnetic 
waves  (e.g.,  laser  beams)  in  dielectric  materials.  This  work  is  a  continuation  of  the  efforts  in 
the  monograph  [2],  where  the  authors  investigated  a  similar  model  using  the  full  Maxwell’s 
equations  together  with  linear  constitutive  relations  in  a  variational  approach  for  the  prop¬ 
agation  of  microwaves  through  dielectric  layers.  The  variational  approach  facilitates  the 
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incorporation  of  antenna  sources  in  the  model  and  provides  an  approximate  computation  of 
the  electromagnetic  transients  like  the  Brillouin  precursor.  We  extended  the  applicability  of 
the  results  in  [2]  in  the  following  sense:  the  model  in  the  monograph  contains  a  linear  consti¬ 
tutive  relationship  describing  material  polarization  in  a  convolution  representation.  Such  a 
formulation  includes  Debye  and  Lorentz  polarization  models  and  is  adequate  to  describe  the 
interaction  of  microwaves  with  dielectric  materials.  However,  the  polarization  mechanism 
cannot  be  assumed  to  be  linear  in  the  case  of  very  high  intensity  electromagnetic  waves, 
like  lasers.  Thus  we  developed  theoretical  results  using  a  representation  of  the  polarization 
by  a  nonlinear  convolution.  In  particular,  we  take  a  physical  model  similar  to  that  in  [2], 
introduce  a  polarization  law 


P(t,z)=f  g(t  —  s,  z)(E(s,  z)  +  f(E(s,  z)))ds,  (2.1) 

and  establish  the  well-posedness  of  the  resulting  system.  This  formulation  can  be  interpreted 
as  a  generalization  of  the  Debye  or  Lorentz  polarization  models  in  the  sense  that  the  polar¬ 
ization  dynamics  is  driven  by  a  nonlinear  function  of  the  electric  field.  We  show  the  global 
existence-uniqueness  and  continuous  dependence  on  data  of  weak  solutions  under  general 
assumptions  on  the  nonlinearity  /.  These  conditions  are  satisfied  by  a  number  of  locally 
polynomial  nonlinearities  proposed  in  the  nonlinear  optics  literature  (e.g.,  see  the  “second 
order”  optical  models  involving  cubic  nonlinearities  in  [10,  11].)  Based  on  these  theoreti¬ 
cal  results,  rigorous  computational  methods  can  be  developed  for  the  forward  problem  and 
subsequently  for  the  associated  inverse  problems.  After  a  detailed  description  of  the  specific 
problem  studied  we  refer  to  [3,  7]  for  the  technical  details  of  the  proofs. 

We  consider  Maxwell’s  equations  applied  to  a  specific  physical  problem  as  depicted  in 
Figure  1.  An  infinite  slab  of  material  is  interrogated  by  a  normally  incident  polarized  plane 
wave  windowed  pulse  originating  at  an  antenna  source  z  =  0  in  free  space  =  [0,  z\].  The 
slab  of  material  in  57  =  [zi ,  z2)  is  assumed  to  be  homogeneous  in  the  directions  orthogonal  to 
the  direction  z  of  propagation  of  the  plane  wave.  Under  these  assumptions  the  strength  of 
the  electric  and  magnetic  fields  in  Q  and  flo  can  be  represented  by  the  scalar  functions  E(t,  z) 
and  H(t,  z),  respectively.  One  can  readily  eliminate  the  magnetic  field  from  the  full  Maxwell 
equations  to  arrive  at  the  strong  formulation  of  the  problem  (for  a  detailed  derivation  see 
[2]) 


Ho tE  +  noIn{z)P  +  fJ'O <?E  -  E"  —  -fio^s  0  <  z  <  z2,  t  >  0, 
IdE  dE\ 


c  dt  dz 


=  0  t  >  0, 

\z= o 

E(t,  z2)  =  0  t  >  0, 

E(0,  z)  =  $(z),  E( 0,  z)  —  \k(z)  0  <  z  <  z2, 


(2.2) 

(2.3) 

(2.4) 

(2.5) 


where  P  denotes  the  polarization  and  Iq  is  the  indicator  function 


In{z) 


0  if  0  <  z  <  zi  1 
1  if  z\  <  z  <  z2,  J  ’ 
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Figure  1:  Geometry  of  the  physical  problem 


where  Z\  <  1  is  the  front  boundary  of  the  slab  and  z2  is  the  (sometimes  unknown  in  in¬ 
terrogation  problems)  back  boundary  (see  Figure  1).  We  note  that  an  absorbing  boundary 
condition  is  placed  at  2  =  0  to  prevent  the  reflection  of  waves.  We  assume  that  there 
is  a  supraconductive  backing  on  the  slab  at  z  =  z2  and  specify  the  boundary  conditions 
accordingly.  Equation  (2.2)  can  be  given  in  the  form 

-E  +  -In(z)P  +  -oE  -  c2E"  =  -- Js> 

€o  Co  Co  Co 

where  e  =  e0(l  +  (er  -  1  )Iq),  and  c2  =  We  denote  ^  by  eT. 

We  assume  that  the  frequency  of  the  interrogating  wave  is  so  high  that  the  dependence 
of  the  polarization  on  the  electric  field  can  no  longer  be  adequately  described  by  a  linear 
constitutive  law.  Specifically,  we  consider  a  polarization  mechanism  that  depends  on  the 
strength  of  the  electric  field  in  the  convolution  in  (2.1).  Thus 

P{t,z)  =  g(t-  s,z)f(E{s,z))ds  +  g{0,z)^f{E{t,z))  +  g(0,z)f{E(t,z)) 

+  [  g(t  -  s,z)E(s,z)ds  +  g(0,z)E(t,z)  +  g(0,z)E(t,z), 

Jo 

where  /  :  1R  — »  M  is  a  nonlinear  (not  necessarily  small)  perturbation  of  the  linear  mechanism. 
This  yields  the  strong  form  of  the  equation 

srE(t,  z)  +  —In(z)(er(z)+g(0,z))E(t,z) 
eo 
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+  —In(z)g(0,z)E(t,z)  +  [  —I{i(z)g(t  —  s,  z)E(s,  z)ds 

eo  eo 

+  -In(z)g(0,z)f(E(t,z))  +  [  —In{z)g{t  -  s,  z)f{E{s,  z))ds 

e0  Jo  e0 

+  -In(z)g(0,z)^f(E(t,z))-c2E"(t,z)  =  --js(t,z),  0  <  z  <  z2.  (2.6) 

eo  at  Co 

In  the  physical  problem  z2  is  assumed  to  be  unknown,  and  it  is  desirable  to  estimate  it  from 
given  data.  Since  the  existence  proof  is  constructive  in  the  sense  that  the  numerical  method 
we  use  to  solve  this  problem  (for  both  forward  and  inverse  problems)  follows  the  theoretical 
arguments,  it  is  desirable  to  convert  the  problem  to  a  fixed  spatial  domain,  e.g.,  [0,1],  as 
explained  in  [2].  We  then  do  not  need  to  update  the  spatial  discretization  and  the  basis 
functions  for  different  z2,  but  the  same  numerical  framework  can  be  used  throughout  the 
optimization  procedure  in  the  inverse  problems. 

To  this  end  we  first  multiply  (2.6)  by  (j),  where  <j>  £  H1  (0,z2)  and  integrate  from  0  to 
z2  using  integration  by  parts  in  the  last  term  in  the  left.  Then  we  introduce  a  change  of 
variable  by  letting 

u  \  -  I  zi  if  0  <  z  <  zi  1 

X~m  ”\  *  +  (*- *i)£=£  J 

i.e.,  h(z)  =  z  +  (C  -  l)(z  -  zi)I[zuZ2](z)  with  C  =  Then  h'(z )  =  1  +  ((-l)IQ(z)  and 
h'(z)  =  1  +  (C  -  lKn(£),  where  (l  =  [zx,  1]  and  we  adopt  the  notation  that  E,  h,  (j>,  g,  etc.,  are 
the  maps  E ,  h,  (j),  g,  etc.,  after  they  have  been  mapped  from  the  domain  [0,  z2]  to  the  domain 
[0, 1].  Using  this  transformation  we  obtain 

zr^E(t,  ■),  <p(z))  +  (~^—yI^(z)(a  +  $(0,  -))E(t,  •),  <p{z)) 
h'(z)  h'(z)  eo 

+  (T^-In(z)§(P>  •)&(!>  ')>  <p(z))  +  (tjtz:  f  ~  ’)&(*>  -)ds>  £(*)> 

h'(z)  eo  h'(z)  Jo  e0 

+  (~^-r-In(z)l(  0,  -)f(E(t,  •)),  (p(z))  +  (J--  [  -  s,  -)f(E(s,  ■ ))ds ,  (p(z)) 

h'(z)  eo  h'(z)  Jo  e0 

+  ■)  ■),#)  +  ck(t,  0)9(0) 

h'(z )  eo  at 

=  (2-7) 

h'(z)  eo 

where  (•,  •)  is  the  L2(0, 1)  inner  product.  This  effectively  maps  our  system  to  a  fixed  reference 
domain  [0, 1]  and  we  use  this  form  to  define  the  weak  solution  of  the  problem.  We  let  H  — 
L2(0, 1),  V  —  HlR( 0, 1)  =  {(f)  e  H1( 0, 1)|^(1)  =  0}  leading  to  the  Gelfand  triple  ([15,  20]) 
V  t->  H  <->-  V*.  We  say  that  E  £  L°°(0,T;  U)  with  E  £  L2(0,T;H),  E  £  L2(0,  T;V*),  is  a 
weak  solution  if  it  satisfies  for  every  ip  €  V 
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(2.8) 


(erE,  p)v-,v  +  (jE,  p)  +  (PE,  p)  +  (f  a(t-  s,  -)E(s,  -)ds,  p) 

JO 

+{Pf{E),  <p)  +  { Jo  a{t  -  s,  -)f(E(s,  ■ ))ds ,  ip)  +  (7  jtf(E),  <p) 
+(c2h'E',  ip')  +  cE(t ,  0)y?(0)  =  (J(t,  •),  V)v\v 


E(Q,z)  =  $(z),  E(0,z)  =  *{z), 


(2.9) 


where  in  (2.8)  we  drop  the  overtilda  on  variables  and  functions  for  simplification  of  notation 
and  for  z  G  [0, 1],  we  define 


a(t,  z) 

m 

7  4) 
7  4) 

J  ( t ,  z) 


1  1 
h'(z )  eo 
1  1 
h'(z)  e 0 
1  1 
h'(z)  Co 
1  1 
h'(i:)  eo 


44)SM)> 
44)5(M)> 
44)44) +  0(04)), 
44)£(0> 5), 


1  1  k, 

J5(t,  z), 


eT  =  ~ 


h'(z)  Co 

1  ~ 
£r. 


h'{z) 

We  note  that  since  >  0  is  a  piecewise  constant  function,  the  inclusion  of  this  term  in 
ck,  /?, 7, 7,  er,  J  is  essentially  equivalent  to  using  a  modified  e0.  In  developing  the  theoretical 
existence  result  we  can  also  make  the  simplifying  assumption  that  er  =  1  without  loss  of 
generality.  (This,  of  course,  is  not  true  for  computational  efforts.)  Our  goal  is  to  establish 
the  existence  and  uniqueness  of  weak  solutions. 

We  make  the  following  assumptions: 

Al)  The  functions  ,0,7,76  L°°(0, 1)  with 

44)1  <  4,  174)1  <  4- 

A2)  The  function  a  is  bounded  on  [0,T]  x  (0,1)  with 

\a{t,z)\  <  4. 
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A3)  g(0,z)  >  0 ,cr(z)  >  0.  (We  note  that  this  assumption  implies  that  7 (2)  >  0,7(2:)  >  0 
for  every  z  G  [0, 1].) 

A4)  The  nonlinear  function  /  :  2R  — >  JR  is  C1,  with  /( 0)  =  0,  and  f'(z)  >  0  for  all  z  G  JR. 
/  is  also  assumed  to  be  affine  at  infinity,  i.e.,  there  exist  constants  R,  LR  and  Kr  such 
that  for  every  |x|  >  R,  \f(x)\  <  LR\x\  +  KR. 

We  proved  the  following  theorem: 

Theorem  2.1  Under  assumptions  Al)-Af)  the  system  (2.8)-(2.9)  has  a  weak  solution  for 
any  $  G  V,  ^  G  H  and  J  E  Hl( 0,T;  V*).  The  weak  solution  depends  continuously  on  initial 
conditions  and  the  source  term,  i.e.,  the  mapping  (<f>,'Ir,  J)  — >  (E,E)  is  continuous  from 
VxHxH\0,T;V*)toC(0,T;H)xL2weak(0,T-H). 

Details  of  the  proof  can  be  found  in  [3]. 

The  theoretical  results  are  an  initial  effort  in  developing  a  foundation  for  nonlinear  elec¬ 
tromagnetics  and  the  needed  computational  methods  for  such.  The  potential  applications  of 
nonlinear  optics  with  light  as  an  information  carrier  are  widespread  and  are  the  basis  of  op¬ 
tical  information  technology  and  the  use  of  optical  fibers.  They  are  of  increasing  importance 
with  the  development  of  materials  (e.g.,  GaAs,  InP,  KNbOs)  which  possess  outstanding  non¬ 
linear  optical  and  electro-optical  properties  [11].  Our  interests  arise  from  the  potential  of 
using  interrogating  input  pulses  in  the  IR  range  (terahertz)  in  imaging  and  more  specifically 
in  detection  algorithms. 

Nonlinear  optical  effects  are  usually  described  through  nonlinear  polarization  laws  [10] 
and  many  of  these  laws  can  be  formulated  in  the  framework  above  where  we  treat  rather 
general  nonlinearities.  These  include  locally  polynomial  nonlinearities  (e.g.,  so-called  third 
order  and  higher  effects  [10,  11])  that  are  bounded  as  the  magnitude  of  the  E  field  saturates 
polarization  mechanisms  (e.g.,  the  material  ’’freezes”  dielectrically  at  high  field  values).  The 
assumption  that  the  nonlinearity  is  increasing  with  the  magnitude  of  the  E  field  is  again  a 
realistic  one. 

We  consider  full  wave  propagation  (as  opposed  to  the  popular  paraxial  approximate 
models  [10]  found  widely  in  the  literature),  even  here  there  are  some  tacit  assumptions 
that  may  result  in  limiting  approximations  to  phenomena  that  actually  occur  in  nonlinear 
materials.  Specifically,  the  one-dimensional  model  formulated  in  this  section  depends  on  the 
tacit  assumption  that  the  polarization  field  P  in  the  dielectric  remains  parallel  to  the  electric 
field  E.  Even  then,  the  usual  Maxwell  equation  V  •  D  =  0  along  with  the  constitutive  law 
D  —  e0 erE  +  fi(P)P  need  not  result  in  V  •  E  —  0.  This  is  important  in  deriving  the  second 
order  form  of  Maxwell’s  equation  where  the  identity  VxVxl  =  V(VT)-  V2E  results  in 
the  simple  Laplacian  only  if  V  ■  E  =  0  or  one  assumes  this  term  is  negligible  as  often  done  in 
nonlinear  optics  ([10]  p.54-60).  Even  with  some  obvious  inadequacies,  the  model  considered 
here  does  facilitate  treatment  of  a  number  of  mechanisms  of  interest,  including  Debye  and 
Lorentz  materials  that  are  driven  by  nonlinear  E  field  inputs,  e.g.,  P  +  ^P  =  f(E). 
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3  Numerical  algorithm  and  simulations 

Our  models  incorporate  different  linear  and  nonlinear  polarization  mechanisms  which  in¬ 
fluence  the  propagation  characteristics.  First  we  study  a  linear  Debye  model  [12],  given 

by 

tP  +  P  =  e0(£s  -  £oo)E,  t>  0,  (3.10) 

and  then  extend  it  to  a  nonlinearly  forced  Debye  and  a  mechanistically  nonlinear  Debye 
model.  Our  computational  framework  is  based  on  a  variational  formulation  of  Maxwell’s 
equations  described  in  the  previous  section. 

Our  focus  here  is  on  understanding  the  propagation  characteristics  of  pulses  with  different 
carrier  frequencies,  and  for  this  purpose  we  set  up  a  simplified  model.  However,  we  note  that 
this  simplified  formulation  can  be  readily  extended  to  treat  interfaces  and  the  more  general 
inverse  problems.  We  consider  an  infinite  (in  the  x  and  y  direction)  slab  of  homogeneous 
material  of  width  i  with  faces  parallel  to  the  xy  plane.  However,  here  we  place  an  antenna 
inside  the  material,  usually  in  the  middle.  The  input  signal  is  a  planar  electromagnetic  wave 
polarized  with  oscillations  in  the  xz  plane  only.  The  electromagnetic  field  E  is  reduced  to 
one  nontrivial  component,  in  the  x  direction,  at  all  points  of  the  slab,  and  it  is  homogeneous 
in  intensity  in  the  x  and  y  directions. 

For  computational  purposes  we  scale  the  time  variable  by  a  factor  of  c  and  polarization 
P  by  a  factor  of  l/e0,  he.,  we  let  E  =  E(ct),P  =  1  /£0P(ct).  We  express  P  from  (3.10) 
and  substitute  it  into  the  weak  formulation  of  the  model.  The  new  equations  in  the  scaled 
variables  are 

{eTE,  <p)  +  ((rja  +  £dX)E,  ip)  +  (A 2P,  ip)  -  (edA2J E,  ip)  +  (E\  <p')  +  y/£^E(0)ip{0) 
+y/FsE{£)p{e)  =  -r){jSl(p)  (3.11) 

(P  +  A  P,  ip)  =  ( £dXE ,  p),  (3.12) 

where  we  introduced  the  notation  £d  —  es  —  £<»,  A  =  T.  ancj  ^  =  ^oC.  We  discretize  the 
problem  in  the  space  variable  using  a  first  order  Galerkin  finite  element  approximation. 
We  divide  the  interval  [0,  i]  into  N  equal  subintervals  at  the  points  z  —  jh,  j  =  0, ...  N, 
where  h  =  £/N,  and  construct  standard  piecewise  linear  splines  (f>f  (z).  The  finite  dimen¬ 
sional  approximating  subspaces  to  V  will  be  taken  to  be  VN  =  |^,  (t>\,  •  •  • ,  •  Now  we 

approximate  E(t,  z )  and  P(t,  z)  in  this  space  as 

E(t,  z)  «  EN(t,  z)  =  J2  ef  (z)>  (3T3) 

j= o 

P(t,z)  «  PN(t,z )  =  fpf  (t)^(z).  (3.14) 

j= o 

We  choose  the  space  of  test  functions  to  be  VN  also,  and  thus  we  find  that  (3.11)-(3.12)  lead 
to 

£rMe  +  (( rja  +  £dX )M  +  y/FsBD)e  +  A 2Mp  +  (K  -  A 2£dM)e  -  -rjJ,  (3.15) 
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Mp  —  —A  Mp  +  EdXMe 


(3.16) 


where  e  =  (e^,  ef, . . . ,  e#)T,  p  =  (Po  ,Pi  ,  ■  ■  ■  ,Pn)T ,  J  =  (Jo,  ■  ■  ■  ,Jn)t ,  with  = 

(js,<f>i),  i  =  0, ...N,  and 

fi 

Mij  =  =  /  (pi(z)(f)j(z)dz , 

J  0 

Kij  =  (4>i4>j)  =  Jo  4l'l(z)(f>'j(z)dz. 

Here  BD  denotes  a  matrix  with  BD00  —  BD^n  =  1,  while  all  other  elements  of  the  matrix 
are  0.  We  now  employ  a  second  order  central  difference  approximation  in  the  time  variable 
to  solve  (3.15).  We  take 


e(tn) 


en+l  _  2en  +  en- 1 

AB 


,  e(tn)  =  en 


Dn+1 


r,n  —  1 


2A  t 


and 


e(tn)  =  en  « 


en+1  +  en  +  e"  1 
4 


Now  the  Galerkin  approximation  (3.15)  reduces  to  a  linear  equation  that  can  be  solved 
for  en+1  given  en,  en~ 1  and  pn.  We  have  initial  condition  e°  =  0  and  we  approximate 
e1  «  —  js(0,  z).  To  obtain  pn  we  solve  (3.16)  using  a  Crank-Nicholson  method.  This 

approximation  method  is  overall  0(h 2)  when  At  =  0(h).  We  also  note  that  the  method  is 
unconditionally  stable. 

We  remark  that  our  model  is  well-posed  theoretically  for  various  different  linear  and  non¬ 
linear  polarization  mechanisms.  The  input  pulse  has  the  form  J4sin3(wt)x[o,t/](^)<^(^  —  zo), 
where  A  is  the  amplitude  of  the  signal,  and  zo  denotes  the  location  where  it  is  generated 
inside  the  material.  We  use  a  few  (typically  4  to  7)  whole  periods  of  this  signal,  i.e.,  tf  is  an 
integer  multiple  of  2tt /ui.  In  the  first  set  of  experiments  we  observed  the  onset  and  propaga¬ 
tion  characteristics  of  the  Brillouin  precursors  in  the  case  when  a  linear  Debye  polarization 
mechanism  is  assumed.  We  find  that  in  the  1-10  GHz  frequency  range  the  peak  amplitude 
of  the  transient  field  is  attenuated  at  a  much  slower  rate  in  relation  to  propagation  distance 
than  the  amplitude  of  the  carrier  frequency  (see  Figures  2  and  3).  Indeed,  a  closer  analysis 
of  the  attenuation  rate  reveals  that  while  the  carrier  frequency  decays  exponentially,  the  at¬ 
tenuation  of  the  transient  is  only  algebraic,  and  proportional  to  approximately  x“0,62  in  the 
1  GHz  case,  and  x-0,59  in  the  10  GHz  case.  This  is  in  accordance  with  the  results  reported 
in  [1]  and  the  theoretical  considerations  based  on  asymptotic  analysis  in  [17]. 

In  the  0.1  to  1  THz  regime  the  carrier  frequency  does  not  propagate  inside  the  material, 
only  the  precursor  enters.  The  attenuation  rate  of  the  amplitude  of  the  leading  transient 
with  respect  to  propagation  length  in  the  1  THz  case  is  approximately  proportional  to  x-1,59. 

Next  we  take  a  nonlinearly  forced  Debye  polarization  model  given  by 


TP  +  P  —  £Q£dE  +  /  ( E ), 


(3.17) 
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Signal  alz*  1.5  (meters) 


t  (ns) 


Signal  at  z*2.1  (meters) 


t(ns) 


Figure  2:  (left)  The  electric  field  recorded  at  the  antenna  and  at  0.6,  0.75  and  1  m  in  the  material 
in  the  first  50  ns  of  propagation,  1  GHz  example  ;  (right)  Peak  of  the  leading  transient  and  the 
carrier  frequency  versus  propagation  length. 


where  f(E)  —  /3E 3.  This  is  an  approximation  to  a  saturation  limited  nonlinearity  required 
by  our  theorem.  We  implemented  the  same  numerical  method  outlined  above,  except  that 
in  this  case  (3.15)-(3.16)  contain  nonlinear  terms.  We  used  a  functional  iteration  in  the 
nonlinear  version  of  (3.15)  to  obtain  e",  which  is  then  used  to  update  the  polarization  termp". 
Comparison  of  the  nonlinearly  forced  model  (3.17)  with  the  corresponding  linear  dynamics 
is  depicted  in  Figures  4  and  5  for  the  10  GHz  and  1  THz  case,  respectively.  In  these 
simulations  we  have  a  weak  nonlinearity  with  ft  =  —  5  x  10~6  and  the  amplitude  of  the  input 
signal  is  small,  A  =  10.  The  linear  and  nonlinear  results  do  not  differ  substantially  in  the  10 
GHz  case  (see  Figure  4).  In  that  case  if  f3  is  positive  we  observed  that  a  nonlinearly  forced 
polarization  dynamics  leads  to  a  signal  whose  main  part  arrives  slightly  earlier  and  is  slightly 
larger  than  the  corresponding  portion  of  the  signal  in  the  linear  material.  More  results  and 
comparisons  can  be  found  in  [7]. 

In  a  set  of  parallel  numerical  experiments  we  implemented  a  nonlinear  Debye  polarization 
dynamics  given  by  rP  +  P  +  sP3  =  e0(£s  —  £<x>)E,  where  s  is  a  small  parameter.  The  cubic 
nonlinear  term  is  chosen  since  most  biological  tissues  are  not  expected  to  have  an  axis 
of  symmetry.  However,  it  must  be  considered  as  an  approximation  to  a  saturation  limited 
nonlinear  mechanism  that  can  be  anticipated  for  orientational  polarization.  Large  departures 
from  the  linear  behavior  can  be  observed  even  for  weak  nonlinearities  if  the  amplitude  of  the 
input  signal  is  sufficiently  large.  In  our  simulations  s  —  10-3  and  the  amplitude  A  of  the 
input  signal  is  in  the  range  1010  to  1012.  Figures  6  and  7  show  the  comparison  between  the 
nonlinear  Debye  model  and  the  corresponding  linear  dynamics  in  the  10  GHz  and  1  THz 
case,  respectively. 

In  the  10GHz  case  we  find  that  initially  the  peak  of  the  transient  is  attenuated  at  a 
faster  rate  than  in  the  linear  case  if  the  amplitude  of  the  input  signal  is  sufficiently  large 
(>  5  •  109).  As  the  propagation  continues,  the  attenuation  rate  seems  to  approximate  that 
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Figure  3:  (left)  The  electric  field  recorded  at  the  antenna  and  at  0.0033,  0.00833  and  0.0166  m 
propagation  in  the  material  in  the  first  2  ns,  10  GHz  example  ;  (right)  Peak  of  the  leading  transient, 
loglog  plot. 
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Figure  4:  Comparison  of  the  linear  (solid  line)  and  nonlinearly  forced  (dashed  line) 
results,  10  GHz  example. 


of  the  linear  case.  Some  of  these  findings  (especially  the  details  of  the  nonlinear  models  and 
corresponding  simulations  and  comparisons  with  the  linear  models)  are  presented  in  [7]. 
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Signal  at  z=0. 0041 6667  (meters) 


Figure  5:  Comparison  of  the  linear  (solid  line)  and  nonlinear ly  forced  (dashed  line) 
results,  1  THz  example. 


4  Models  for  polymeric  materials 

As  one  component  in  the  investigation  of  the  interaction  of  electromagnetic  waves  with 
viscoelastic  materials  we  further  developed  our  models  for  large  dynamic  deformations  of  a 
polymeric  rod.  These  materials  exhibit  large  hysteresis  loops  between  the  stress  and  strain, 
which  is  usually  captured  by  nonlinear  constitutive  relationships  and  is  thought  of  as  the 
manifestation  of  a  ‘memory’  in  the  material.  In  contrast  to  our  earlier  efforts  based  on  a 
pseudo-phenomenological  description  of  this  phenomenon  [4,  9]  we  developed  a  molecular, 
i.e.,  physics  based,  model  that  can  incorporate  the  multiscale  aspect  of  the  problem.  This 
model  is  based  on  a  nonlinear  extension  of  the  linear  reptation  models  of  Doi  and  Edwards 
[13]  and  the  work  of  Johnson  and  Stacer  [14].  We  assume  that  the  material  consists  of 
physically-constrained  and  chemically  cross-linked  molecules  whose  motion  with  respect  to 
each  other  constitutes  an  ‘internal  dynamics’  in  the  model  that  contributes  to  the  overall 
system  response.  We  showed  in  [5]  and  [6]  that  the  earlier  models  of  [4,  9]  are  special 
cases  of  this  formulation  if  one  considers  uniform  (or  a  finite  number  of)  internal  dynamics. 
In  the  new  model  we  treat  a  continuum  of  internal  dynamics  that  are  averaged  according 
to  a  probability  distribution.  The  resulting  system  is  a  probability  measure  dependent 
partial  differential  equation  for  which  we  established  well-posedness  in  [8].  This  approach 
treats  hysteresis  as  a  multiscale  phenomenon  and  is  shown  to  lead  to  computationally  useful 
approximations  [8]. 
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Figure  6:  Comparison  of  the  linear  (solid  line)  and  nonlinear  (dashed  line)  results, 
10  GHz  example. 
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